/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		tensor.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-02-18 13:58
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:23
*
*   Description:
*
=============================================================================*/
#if !defined(BASE_TENSOR_H)
#define BASE_TENSOR_H

#include "../geometry/shape/dimension_type.hpp"
#include "basic_arithmetic.hpp"
#include <assert.h>
#include <math.h>

#include <iostream>

namespace geoxfem {
namespace arithmetic {

struct stress {};
struct strain {};

template <typename DimensionType, typename tensor_type>
struct Tensor : Matrix3R {

  /*constructor{{{*/
  Tensor() : Matrix3R{Matrix3R::Zero()} {}

  template <typename OtherDerived>
  Tensor(const Eigen::MatrixBase<OtherDerived>& other) : Matrix3R(other) {
    assert(other.cols() == 3);
    assert(other.rows() == 3);
  }
  /** This method allows you to assign Eigen expressions to MyVectorType
   *
   * @param other
   * @return
   */
  template <typename OtherDerived>
  Tensor& operator=(const Eigen::MatrixBase<OtherDerived>& other) {
    Matrix3R::operator=(other);
    return *this;
  }
  /*}}}*/
  /*convert to voigt{{{*/
  /** init from voigt
   *
   * @param voigt  :
   * @param _type
   */
  void fromVoigt(const VectorXR&) {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  /** return Tensor form
   *
   * @param _type
   * @return
   */
  VectorXR Voigt() const {}
  /*}}}*/
  /*deviatorics{{{*/
  _real dm() const {
    return this->trace() / _real(3);
  }
  _real s11() const {
    return this->operator()(0, 0) - dm();
  }
  _real s22() const {
    return this->operator()(1, 1) - dm();
  }
  _real s33() const {
    return this->operator()(2, 2) - dm();
  }
  /*}}}*/
  /*invariants and their derivatives{{{*/

  /** get invariant 1. I1=d11+d22+d33
   *
   * @return
   */
  _real I1() const {
    return this->trace();
  }
  /** get invariant 2. I2=d11*d22+d22*d33+d11*d33-d12^2-d13^2-d23^2
   *
   * @return
   */
  _real I2() const {
    return this->operator()(0, 0) * this->operator()(1, 1)
           + this->operator()(1, 1) * this->operator()(2, 2)
           + this->operator()(0, 0) * this->operator()(2, 2)
           - this->operator()(0, 1) * this->operator()(0, 1)
           - this->operator()(1, 2) * this->operator()(1, 2)
           - this->operator()(0, 2) * this->operator()(0, 2);
  }
  /** get invariant 3. I3= d11*d22*d33 + 2d12*d13*d23 - d11*d23*d23 -
   * d22*d13*d13 - d33*d12*d12
   *
   * @return
   */
  _real I3() const {
    return this->determinant();
  }
  /** get the derivatives of I1 with respect to the component
   *
   * @return
   */
  VectorXR dI1dS() const {
    VectorXR dI1dS(6);
    dI1dS << 1, 1, 1, 0, 0, 0;

    return dI1dS;
  }
  /** get the derivatives of I2 with respect to the component
   *
   * @return
   */
  VectorXR dI2dS() const {
    VectorXR dI2dS(6);

    dI2dS << this->operator()(1, 1) + this->operator()(2, 2),
      this->operator()(0, 0) + this->operator()(2, 2),
      this->operator()(0, 0) + this->operator()(1, 1),
      -2 * this->operator()(0, 1), -2 * this->operator()(0, 2),
      -2 * this->operator()(1, 2);

    return dI2dS;
  }
  /** get the derivatives of I3 with respect to the component
   *
   * @return
   */
  VectorXR dI3dS() const {
    VectorXR dI3dS(6);

    dI3dS << this->operator()(1, 1) * this->operator()(2, 2)
               - this->operator()(1, 2) * this->operator()(1, 2),
      this->operator()(0, 0) * this->operator()(2, 2)
        - this->operator()(0, 2) * this->operator()(0, 2),
      this->operator()(0, 0) * this->operator()(1, 1)
        - this->operator()(0, 1) * this->operator()(0, 1),
      2 * this->operator()(0, 2) * this->operator()(1, 2)
        - 2 * this->operator()(0, 1) * this->operator()(2, 2),
      2 * this->operator()(0, 1) * this->operator()(1, 2)
        - 2 * this->operator()(0, 2) * this->operator()(1, 1),
      2 * this->operator()(0, 1) * this->operator()(0, 2)
        - 2 * this->operator()(1, 2) * this->operator()(0, 0);

    return dI3dS;
  }
  /*}}}*/
  /*deviatoric invariant and their derivatives{{{*/
  /** get deviatoric tensorpacketp
   *
   * @return
   */
  const Tensor getDeviatoric() const {
    Tensor Deviatoric = *this;
    Deviatoric -= this->Identity() * this->trace() / 3.;

    return Deviatoric;
  }
  /** get deviatoric tensor invariant 1
   *
   * @return
   */
  _real J1() const {
    return _real(0);
  }
  /** get deviatoric tensor invariant 2
   *
   * @return
   */
  _real J2() const {
    /// also write as
    /// 1/6*[(d11-d22)^2+(d22-d33)^2+(d11-d33)^2]+d12^2+d13^2+d23^2
    const _real _I1 = I1(), _I2 = I2();
    return _I1 * _I1 / _real(3) - _I2;
  }
  /** get deviatoric tensor invariant 3
   *
   * @return
   */
  _real J3() const {
    const _real _I1 = I1(), _I2 = I2(), _I3 = I3();
    return _real(2) * _I1 * _I1 * _I1 / _real(27) - _I1 * _I2 / _real(3) + _I3;
  }
  /** get the derivatives of p with respect to the component
   *
   * @return
   */
  MatrixXR dpdS() const {
    return dI1dS() / 3.;
  }
  /** get the derivatives of J2 with respect to the component
   *
   * @return
   */
  VectorXR dJ2dS() const {
    VectorXR dJ2dS(6);

    dJ2dS << s11(), s22(), s33(), 2 * this->operator()(0, 1),
      2 * this->operator()(0, 2), 2 * this->operator()(1, 2);

    return dJ2dS;
  }
  /** get the derivatives of J3 with respect to the component
   *
   * @return
   */
  VectorXR dJ3dS() const {
    Matrix3R s = getDeviatoric();
    Matrix3R _dJ3dS;
    VectorXR dJ3dS(6);

    _dJ3dS = s * s - Matrix3R::Identity() * J2() * _real(2) / _real(3);
    dJ3dS << _dJ3dS.operator()(0, 0), _dJ3dS.operator()(1, 1),
      _dJ3dS.operator()(2, 2),
      _dJ3dS.operator()(0, 1) + _dJ3dS.operator()(1, 0),
      _dJ3dS.operator()(0, 2) + _dJ3dS.operator()(2, 0),
      _dJ3dS.operator()(1, 2) + _dJ3dS.operator()(2, 1);

    return dJ3dS;
  }
  /*}}}*/
  /*misc{{{*/

  /** get /sigma_{Octahedral}
   *
   * @param /  :
   * @return
   */
  _real Soct() const {
    return std::sqrt(_real(2) / _real(3) * J2());
  }
  /** get mean stress
   *
   * @return
   */
  _real p() const {
    return this->trace() / _real(3);
  }
  /** get deviatoric stress
   *
   * @param J2() :
   * @return
   */
  _real q() const {
    return std::sqrt(_real(3) * J2());
  }
  /** get Von Mises stress
   *
   * @param J2() :
   * @return
   */
  _real Mises() const {
    return std::sqrt(_real(3) * J2());
  }
  /**
   * @brief effective stress
   *
   * @return
   */
  _real Sbar() const {
    return std::sqrt(_real(3) * J2());
  }
  /** get the derivatives of q with respect to the component
   *
   * @return
   */
  VectorXR dqdS() const {
    if(J2() == 0)
      return VectorXR::Zero(6);
    return std::sqrt(_real(3) / J2()) * dJ2dS() / _real(2);
  }
  /** get the derivatives of Mises with respect to the component
   *
   * @return
   */
  VectorXR dMisesdS() const {
    if(J2() == 0)
      return VectorXR::Zero(6);
    return std::sqrt(_real(3) / J2()) * dJ2dS() / _real(2);
  }
  /** get the derivatives of rho with respect to the component
   *
   * @return
   */
  VectorXR dRhodS() const {
    if(J2() == 0)
      return VectorXR::Zero(6);
    return std::sqrt(_real(2) / J2()) * dJ2dS() / _real(2);
  }
  /*}}}*/
  /*Haigh-Westergaard coordinates{{{*/

  /** get lode theta. anti-clockwise. start from x-axis direction at PI plane
   * domain: (-PI/6,PI/6)
   *
   * @return
   */
  _real lode2() const {
    if(J2() == 0)
      return 0;
    return _real(1) / _real(3)
           * std::asin(-std::pow(_real(3) / J2(), _real(1.5)) * J3()
                       * _real(0.5));
  }
  /** get  Haigh-Westergaard coordinates Xi
   *
   * @return
   */
  _real Xi() const {
    return I1() * root3Reciprocal;
  }
  /** get  Haigh-Westergaard coordinates rho
   *
   * @param J2() :
   * @return
   */
  _real rho() const {
    return std::sqrt(_real(2) * J2());
  }
  /**  get Haigh-Westergaard coordinates lode angle. clockwise. start from 150
   * degree direction at PI plane domain: (-PI/6,PI/6)
   *
   * @return
   */
  _real lode() const {
    if(J2() == 0)
      return 0;
    return _real(1) / _real(3)
           * std::acos(std::pow(_real(3) / J2(), _real(1.5)) * J3() / _real(2));
  }
  /** get the derivatives of lode lode with respect to sigma
   *
   * @return
   */
  VectorXR dLodedS() const {
    const _real _lode = lode();
    const _real _J2   = J2();
    const _real _J3   = J3();

    return root3 / _real(2) / std::sin(_real(3) * _lode)
           / std::pow(_J2, _real(1.5))
           * (_real(1.5) * _J3 / _J2 * dJ2dS() - dJ3dS());
  }
  /*}}}*/
  /*principal stress{{{*/

  /** get maximum principal stress
   *
   * @return
   */
  _real Pmax() const {
    return Xi() * root3Reciprocal
           + root2 * root3Reciprocal * rho() * std::cos(lode());
  }
  /** get mid principal stress
   *
   * @return
   */
  _real Pmid() const {
    return Xi() * root3Reciprocal
           + root2 * root3Reciprocal * rho()
               * std::cos(lode() - _real(2) / _real(3) * _real(M_PI));
  }
  /** get minimum principal stress
   *
   * @return
   */
  _real Pmin() const {
    return Xi() * root3Reciprocal
           + root2 * root3Reciprocal * rho()
               * std::cos(lode() + _real(2) / _real(3) * _real(M_PI));
  }
  /*}}}*/
};

/*convert from voigt vector{{{*/
template <>
void Tensor<geometry::plane_stress, stress>::fromVoigt(const VectorXR& voigt) {
  assert(voigt.size() == 3);
  *this << voigt[0], voigt[2], 0,  //
    voigt[2], voigt[1], 0,         //
    0, 0, 0;
}
template <>
void Tensor<geometry::plane_strain, stress>::fromVoigt(const VectorXR& voigt) {
  assert(voigt.size() == 3);
  *this << voigt[0], voigt[2], 0,  //
    voigt[2], voigt[1], 0,         //
    0, 0, 0;
}
template <>
void Tensor<geometry::three_dimension, stress>::fromVoigt(
  const VectorXR& voigt) {
  assert(voigt.size() == 6);
  *this << voigt[0], voigt[3], voigt[4],  //
    voigt[3], voigt[1], voigt[5],         //
    voigt[4], voigt[5], voigt[2];
}
template <>
void Tensor<geometry::plane_stress, strain>::fromVoigt(const VectorXR& voigt) {
  assert(voigt.size() == 3);
  *this << voigt[0], voigt[2] * r_5, 0,  //
    voigt[2] * r_5, voigt[1], 0,         //
    0, 0, 0;
}
template <>
void Tensor<geometry::plane_strain, strain>::fromVoigt(const VectorXR& voigt) {
  assert(voigt.size() == 3);
  *this << voigt[0], voigt[2] * r_5, 0,  //
    voigt[2] * r_5, voigt[1], 0,         //
    0, 0, 0;
}
template <>
void Tensor<geometry::three_dimension, strain>::fromVoigt(
  const VectorXR& voigt) {
  assert(voigt.size() == 6);
  *this << voigt[0], voigt[3] * r_5, voigt[4] * r_5,  //
    voigt[3] * r_5, voigt[1], voigt[5] * r_5,         //
    voigt[4] * r_5, voigt[5] * r_5, voigt[2];
}
/*}}}*/
/*convert to voigt vector{{{*/
template <>
VectorXR Tensor<geometry::plane_stress, stress>::Voigt() const {
  VectorXR       voigt(3);
  voigt << this->operator()(0, 0),  //
    this->       operator()(1, 1),  //
    this->       operator()(0, 1);
  return voigt;
}
template <>
VectorXR Tensor<geometry::plane_strain, stress>::Voigt() const {
  VectorXR       voigt(3);
  voigt << this->operator()(0, 0),  //
    this->       operator()(1, 1),  //
    this->       operator()(0, 1);
  /* this->       operator()(2, 2); */
  return voigt;
}
template <>
VectorXR Tensor<geometry::three_dimension, stress>::Voigt() const {
  VectorXR       voigt(6);
  voigt << this->operator()(0, 0),  //
    this->       operator()(1, 1),  //
    this->       operator()(2, 2),  //
    this->       operator()(0, 1),  //
    this->       operator()(0, 2),  //
    this->       operator()(1, 2);
  return voigt;
}
template <>
VectorXR Tensor<geometry::plane_stress, strain>::Voigt() const {
  VectorXR       voigt(3);
  voigt << this->operator()(0, 0),  //
    this->       operator()(1, 1),  //
    this->operator()(0, 1) + this->operator()(1, 0);
  return voigt;
}
template <>
VectorXR Tensor<geometry::plane_strain, strain>::Voigt() const {
  VectorXR       voigt(3);
  voigt << this->operator()(0, 0),  //
    this->       operator()(1, 1),  //
    this->operator()(0, 1) + this->operator()(1, 0);
  /* this->                         operator()(2, 2); */
  return voigt;
}
template <>
VectorXR Tensor<geometry::three_dimension, strain>::Voigt() const {
  VectorXR       voigt(6);
  voigt << this->operator()(0, 0),              //
    this->       operator()(1, 1),              //
    this->       operator()(2, 2),              //
    this->operator()(0, 1) + operator()(1, 0),  //
    this->operator()(0, 2) + operator()(2, 0),  //
    this->operator()(1, 2) + operator()(2, 1);
  return voigt;
}
/*}}}*/

template <typename DimensionType>
using Stress = Tensor<DimensionType, stress>;
template <typename DimensionType>
using Strain = Tensor<DimensionType, strain>;

}  // namespace arithmetic

template <typename dimensionType, typename tensor_type>
struct traits::dimension_type<
  arithmetic::Tensor<dimensionType, tensor_type>> { /*{{{*/
  using type = dimensionType;
}; /*}}}*/

}  // namespace geoxfem
#endif
